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We investigate noise-induced pattern formation in a model of cancer 
growth based on Michaelis-Menten kinetics, subject to additive and multi- 
plicative noises. We analyse stability properties of the system and discuss 
the role of diffusion and noises in the system's dynamics. We find that ran- 
dom dichotomous fluctuations in the immune response intensity along with 
Gaussian environmental noise lead to emergence of a spatial pattern of two 
phases, in which cancer cells, or, respectively, immune cells predominate. 

PACS numbers: 05.40.-a, 87.10.+e, 89.75.Kd 

1. Introduction 

The study of population dynamics covers a wide range of fields such as 
ecology, cellular and molecular biology and medicine (see e.g. ^ |U E3 E] ) • 
The models of population growth based on nonlinear ordinary differential 
equations, despite their simplicity, can often capture the essence of complex 
biological interactions and explain characteristics of proliferation phenom- 
ena. However, biological processes are not purely deterministic: systems 
existing in nature are subject to natural noises. 

The presence of noise in biological systems gives rise to a rich variety 
of dynamical effects. Random fluctuations may be regarded not only as 
a mere source of disorder but also as a factor which introduces positive 
and organizing rather than disruptive changes in the system's dynamics. 
Some of the more important examples of noise induced effects are stochastic 
resonance [SJ |2] , resonant activation H3 El , noise-induced transitions |3 
ITU] , noise-enhanced stability and pattern formation [Tl 01 IT2]. 

The effect of cell-mediated immune surveillance against cancer [131 1141 
ITS] may be a specific illustration of the coupling between noise and a bio- 
logical system. Most of tumoral cells bear antigens which are recognised as 
strange by the immune system. A response against these antigens may be 
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mediated either by immune cells such as T-lymphocytes or other cells, not 
directly related to the immune system (like macrophages or natural killer 
cells) . The process of damage to tumour proceeds via infiltration of the lat- 
ter by the specialised cells which subsequently develop a cytotoxic activity 
against the cancer cell-population. The series of reactions between the cyto- 
toxic cells and the tumour tissue may be considered to be well approximated 
by a saturating, enzymatic-like process whose time evolution equations are 
similar to the standard Michaelis-Menten kinetics |131 116j . Random vari- 
ability of kinetic parameters defining that process may effect the extinction 
of the tumour. 

In our previous article [B] we discussed the noise-induced effect of res- 
onant activation in a spatially homogeneous model of cancer growth 
based upon the above-mentioned kinetic scheme. In the present paper, we 
focus on the study of a spatially inhomogeneous system, namely, we inves- 
tigate how global environmental noise as well as fluctuations in the immune 
response parameter effect the emergence of spatial patterns. 

2. The Model 

The interaction between cancer cells and cytotoxic cells will be described 
by use of the predator-prey model based upon the Michaelis-Menten kinetic 
scheme [HI H3 ECU E3 El EE] This model is a classical one and has been exten- 
sively studied since the 1970s. Its validity has been verified experimentally 
e.g. in |19j . where the authors examined the mechanism of immune rejection 
of a tumour induced by Moloney murine sarcoma virus. The behaviour of 
the cellular populations may be represented by the following scheme: 

First, the cytotoxic cells bind to the tumour cells at rate k\\ subse- 
quently, the cancer cells which have been bound are killed and the complex 
dissociates at rate ^2; finally, dead cancer cells decay at a rate k%. The 
process can be described schematically: 

X + Y^Z^Y + P-^Y. (1) 

X represents here the population of tumour cells. Y, Z and P are popu- 
lations of active cytotoxic cells, bound cells and dead tumour cells, respec- 
tively. Following the original presentation ^U], we assume that (i) cancer 
cells undergo replication at a rate proportional to a time constant A; (ii) as 
a result of cellular replication in limited volume, a diffusive propagation of 
cancer cells is possible, with a transport coefficient proportional to the repli- 
cation rate and local density of cancer cells; (iii) dead cancer cells undergo 
elimination at rate £3; (iv) free cytotoxic cells can move with a "diffusion" 
coefficient D. The spatio-temporal evolution of the tumour due to the above 
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processes can be then described by a set of balance equations: 

' J = A[l - ( x + p)]x- k iy x + A(l -p)V 2 x + XxV 2 p 
§ = -k lV x + k 2 z + DX7 2 y ^ 
P = fciyx - k 2 z 
-Jl = k 2 z - k 3 p 

The x(~r*, t), y(!*, t), z{j*, t) and p(T*, t) are local densities of cells at point 
"r\ Finally, we impose an additional condition on the model: the total 
number of active and bound cytotoxic cells should remain constant: 

y(t) + z(t) = const. = E. (3) 

Below, we analyse several versions of the model: 

• Bifurcation analysis, based on a bare kinetics model without spatial 
diffusion (Sec. ETTj) 

• Incorporation of Fickian diffusion terms in evolution equations for 
x(T > , t), p(/r*, t) respectively, which leads to a wavefront solution (Sec. 
El 

• We determine the effect of the dichotomous switching in the kinetic 
parameter k\ and discuss its role also for the model system in which 
each of the kinetic equations is additionally driven by a Gaussian noise 
term of the same intensity (Sec. 14.31 14.41 14. 5|) 

• At the last step of the model complexity, we analyse how the joint 
effect of diffusion, additive Gaussian fluctuations and independent di- 
chotomic switching in k\ parameter affect the cancer cells population 
dynamics (Sec. I4.6|) . 

3. Simulation 

The basic aim of this work was to study the behaviour of the system © 
subject to additive and multiplicative noises: 

' § = \[l-{x+p))x-(k l +r,(t))yx+ 

+A(1 - p)V 2 x + XxV 2 p + <r£(y, t) 
< |f = -(h +r ] {t))yx + k 2 z + DV 2 y + o£{r,t) (4) 

p = {h + 7](t))yx - k 2 z + v£C?,t) 
k = k 2 z-k 3 p + a£(~r>,t) 

The multiplicative dichotomous Markovian noise n(t) = ±A with mean 
frequency 7 and autocorrelation < rj{t)rj{t') >= -^e -27 '* - *'' models fluc- 
tuations in immune response. The additive Gaussian noise £,(^,t) with 
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autocorellation < £(!*, t)£(~r* , t') >= 5(1* — 1* )6(t — t') depicts external 
environmental fluctuations. We assume that its intensity a is same for each 
variable of the system. Both noises are assumed to be statistically indepen- 
dent: < £(~f*,i)7/(s) >= 0. 

3.1. Numerics 

We have solved the stochastic differential equations © numerically, us- 
ing the Euler scheme, on a 128 x 128 square lattice. According to the 
statistical properties of T](t), the waiting time between two switchings was 
generated from the exponential distribution. 

Since x, y, z and p are densities, their values never can be greater than 1 
or less than 0. Consequently, following boundary conditions have been im- 
posed on the simulated system: if x, y, z or p becomes less than or greater 
than 1 at a given time step, we adjust its value to or to 1, respectively. 

3.2. Simulation Results 
We performed a simulation with the following values of parameters: 

A = 0.5, D = 0.05, a = 0.01, A = 0.5, , , 

h = 1.75, k 2 = 0.1, k 3 = 0.1, 7 = 0.01. ^ 

7 is the mean rate of switching in rj(t). Initial conditions: 

x(T,0) = 0, 0) = 0.4, z(T,0) = 0, p( r\0) = 0. (6) 

The values of parameters and initial conditions have been chosen so that 
we could obtain a distinct pattern: The immune response rates k\ + A and 
k± — A, along with A lead to two different types of stationary behaviour 
(see Sec. 14.5(1 . The mean switching rate is one or two orders of magnitude 
slower than other kinetic parameters, which gives the system a possibility 
to approach the stationary states. The environmental noise intensity a has 
been chosen in such a way that it allows the system to jump between both 
mentioned states. The noise is, however, weak enough to let the system form 
a pattern. At the selected value of D, the pattern has sufficiently distinct 
boundaries and is relatively stable (at higher values of D it would dissolve 
quickly, whereas at smaller D it would form small "grains"). Parameters k 2 
and &3 have been chosen by a trial-and-error procedure: k% is responsible 
for the dissociation of z into y + p. If the dissociation rate is large, then 
the active immune cells y are being released faster and thus the immune 
response is more effective. The &3 parameter determines the rate at which 
dead cancer cells are eliminated. Since dead cells occupy the living space, 
this parameter controls the effective replication rate of cancer cells. 
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Fig. 1. Snapshots from the simulation of the system (@J. Time evolution of the 
spatial distribution of x, y and z. Time counter in upper left corner of each image. 
Simulation time T — 10000. Initial conditions: x = 0, y = 0.4, z — 0,p — every- 
where. Parameters: A = 0.5, D = 0.05, a = 0.01, A = 0.5, ki = 1.75, k 2 = 0.1, k 3 = 
0.1,7 = 0.01. Light pixels: high concentration, dark pixels: low concentration of 
cells. 



The simulation results are presented in Fig. ^ After some time, we 
observe the emergence of the "y-phase" islands (where the immune cells y 
predominate) within the relatively homogeneous "x-phase" (in this phase 
cancer cells x prevail). The phase boundaries move back and forth depend- 
ing on the dichotomous changes in the immune response intensity k\ ± A. 



4. Analysis 

In order to explain the behaviour of the simulated system, we will analyse 
its stability properties as well as the role of diffusion, additive noise and 
dichotomous noise in the k\ parameter. 
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4-1. Stability analysis 
The stationary points {x* ,y* , z* ,p*} of the system are given by: 

{0,y*, 0,0} (7) 

and 

f fc 3 A - hy* + hk 3 y* A - hy* k x y* A - kiy* \ 

I A k 3 + k iy * ,V ' AA: 2 h + hy*' A h + hy*)' U 

According to El the value of y* is arbitrary. However, if the conditional 
is involved, it defines the value of y*: 

y* + z* = const. = E. (9) 

The sets of stationary points form two branches in the x-y-z-p-space. 
The branch changes its stability at the point {0, ^-,0,0}. It is repelling 
for < y < and attracting for ^- < y < 1. The branch (jSJ is attracting 

for0<y<^(-l + yiT^). For^(-l + v /TT^) < y < 1 it consists 
of saddle points. 

Trajectories of the deterministic system can lie only on the plane @, 
given by the initial conditions for y and z. Hence, the stationary points 
of such a system lie on the intersection of the plane @ and the branches 
0, ©■ Depending on the values of parameters and initial conditions, the 
system can have 1 (attracting), 2 (attracting and repelling) or 3 (attracting, 
saddle and attracting) stationary points (see Fig. [21 [SJ)- 

When E, k2, k 3 are fixed, the stability properties of the system depend 
on the ki parameter (see Fig. |1J, i.e. the immune response efficiency, which, 
in our simulation, was controlled by the dichotomous noise. 

In the next subsections, we will analyse how the presence of noises and 
diffusion affects the behaviour of the system. 

4-2. Spatially inhomogeneous system without noise 

The solutions to the deterministic system with diffusive terms have the 
form of travelling wavefronts j5U] . 

4-3. Spatially homogeneous system with dichotomous multiplicative noise 

After the introduction of the dichotomous noise rj(t) = ±A into the k\ 
parameter, trajectories of the system lie on the planes y + z = E and wander 
between two stationary points according to the current value of k% ± A (see 
Fig. El). 
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Fig. 2. Sets of stationary points JJJ, (JHJ) and their stability (shown only schemat- 
ically by arrows), x-y-z-projection. Dashed lines: planar projections of the curve. 
Parameters: A = 0.5, ki = 1.25, k 2 = 0.1, fc3 = 0.1 




o 



Fig. 3. Sets of stationary points 07© m y-z-p-projection. Dashed lines: planar 
projections of the curve. Parameters: A = 0.5, k\ — 1.25, k 2 — 0.1, k3 = 0.1 
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Fig. 4. Stationary points (open circles) lie on the intersection of the plane y + z = E 
and the curves given by {7J|,|(8J|. As k\ grows, the number of intersections changes 
from 2 to 3 and, subsequently, from 3 to 1. (An arrow shows schematically the 
direction in which the intersection points move as k\ grows.) The values of the 
remaining parameters are here: E = 0.6, A = 0.6, = 0.1, k% = 0.1 




Fig. 5. Example trajectories of the system with dichotomous noise only, for two 
different values of E = y + z: E = 0.9 and E — 0.2. Parameters: A = 0.5,7 — 
0.05, ki = 0.5, A = 0.25, k 2 = 1, k 3 = 1. Simulation time T = 700. 
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Fig. 6. Two example trajectories of the spatially homogeneous system with multi- 
plicative dichotomous noise, with additive Gaussian noise, starting from two dif- 
ferent initial conditions: y + z = E = 0.9 and y + z = E = 0.2, in x-y-projection. 
(7 = 0.002, other parameters same as in Fig. [S] To compare with the trajectories 
without additive noise, see Fig. |3] 

4-4- Spatially homogeneous system with dichotomous multiplicative noise 



In the next step of complexity, we add the term cr^T^i) (Gaussian 
white noise) to each equation, assuming that the additive noise acts in the 
same way on each variable of the system. Here, the trajectories do not stay 
on constant planes y + z = E any more. Due to the additive noise, they 
"slip off" from their initial planes (see FigH3[7|). 



Let us examine a system with the additive Gaussian noise and diffusion, 
but without the dichotomous noise. Instead, we will take into consideration 
two situations where the intensity of the immune response is constant: 



and additive Gaussian white noise 



4-5. System with Gaussian noise and diffusion 



Ki = jfei - A 



(10) 



and 



K 1 = k 1 + A. 




The value of A is same as in the first simulation (Sec. I3.2|) . i.e. it corresponds 
to one of the dichotomous noise states. 
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For K\ = k\ — A, the "x-phase" is more stable. The trajectory starts 
with initial conditions © which is exactly the point {0, ^-,0,0} where two 
branches of stationary points cross. It is then very likely that the trajectory 
falls down onto the lower branch of the attractor and stays in its neighbour- 
hood (see Fig. I3ITU1) . 

For K\ = ki + A, the "y-phase" is more stable. Starting far away from 
the lower branch of the attractor, the trajectory remains close to the upper 
branch. Moreover, it climbs higher and higher because of the reflecting 
boundary at x = 0, y = 0, z = 0,p = 0, x = l,y = 1, z = l,p = 1 (see 
Fig. 1111121) . This condition imposed on the boundaries is justified by the 
requirement that the population density cannot be negative nor greater than 
1. 

One can observe a synchronisation effect: the values of y, z and p de- 
crease when x increase, and vice versa (see Fig. HUI12j) . One can also 
notice the stabilising effect of diffusion. Fig. EJcompares trajectories of two 
systems: with and without diffusion, driven by additive Gaussian noise. 
Without diffusion, the trajectory wanders up and down along the branch of 
stationary points. The trajectory stabilised by diffusion stays longer in the 
neighbourhood of its initial plane E. 

4-6. System with Gaussian noise, dichotomous noise and diffusion 

The system @ is a combination of all cases analysed above. Dichoto- 
mous noise switches the system between states where either the "x-phase" 
or the "y-phase" is preferred. This causes the emergence of separate "is- 
lands" of these two phases. Their boundaries move due to diffusion and the 
direction and speed of the motion depends on the current value of k\ + rj(t) 
(see Fig. I1I15|) . In the x-y-z-space, the trajectories climb up towards the 
region where y is close to 1 (the upper part of the y-phase attractor). This 
"climbing" effect is caused by the boundary conditions imposed on the sys- 
tem: Since the same positive or negative value of cr£(~r*, t) is being added 
to each equation, there is a certain preferred direction in which the the 
trajectory moves. It cannot, however, cross the boundaries and thus, near 
the attracting branch of ((7J), the average direction of motion is "upwards" 
(i.e. towards the increasing values of y) because the motion in the opposite 
direction is blocked due to the boundary conditions (see Fig |14|) . 

When trajectories get to the upper part of the "y-phase" attractor, they 
remain in its neighbourhood for all time because the "x-phase" attractor is 
too distant. This distance is determined by the choice of parameters, namely 
by the position of point {0, ^-,0,0} where the other branch of stationary 
points begins (see Fig. I14j) . This effect is the reason why the "y-phase" 
finally spreads all over accessible space in our simulation. 
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Fig. 7. Two example trajectories of the spatially homogeneous system with multi- 
plicative dichotomous noise and additive Gaussian noise, starting with two different 
initial conditions: E = 0.9 and E — 0.2, in y-z-projection. The trajectories do not 
stay on their initial planes E (denoted by sloped lines), but move up or down due 
to the Gaussian noise. Curves: branches of stationary points for switching value 
of k\ ± A . a = 0.002, other parameters same as in Fig. |SJ 
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Fig. 8. Stabilising effect of diffusion. Gray: trajectory of a system with diffusion 
(D = 0.5) and additive Gaussian noise, recorded at point [20, 20] on the spatial 
lattice. Black: trajectory of a system with additive Gaussian noise, but without 
diffusion (D = 0). Parameters: A = 0.5, a = 0.005, k x = 0.25, k 2 = 1, k 3 = 1. 
Simulation time T = 3000. 
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Fig. 9. An example trajectory of the system with additive Gaussian noise and 
diffusion, recorded at point [20, 20] on the spatial lattice. The immune response 
parameter K\ = k\ — A = 1.25. The trajectory falls down onto the "x-phase" 
attractor. A = 0.5, D = 0.05, a = 0.01, A = 0, k 2 = 0.1, k 3 = 0.1. Initial conditions: 
x = 0, y = 0.4, z = 0,p = 0. Straight lines: the profiles of example E planes. 

5. Conclusions 

We have performed a simulation of a spatially inhomogeneous model of 
cancer growth with additive Gaussian noise and multiplicative dichotomous 
noise. The multiplicative noise controls the efficiency of the immune re- 
sponse whereas the external environmental fluctuations have been modelled 
by the additive Gaussian noise. 

The presence of noise in biological systems may be regarded not only as 
a mere source of disorder but also as a factor which introduces positive and 
organising rather than disruptive changes in the system's dynamics: In our 
model, we find that the presence of a global (i.e. depending only on time) 
multiplicative dichotomous noise in a system perturbed by a spatially inho- 
mogeneous additive Gaussian noise leads to emergence of a spatial pattern 
of two "phases". The "phases" are distinct areas in which cancer cells, or, 
respectively, immune cells predominate. The pattern is not stable: domain 
boundaries move due to a diffusion effect, and the direction (and speed) of 
that motion is determined by the current value of the dichotomous noise, 
i.e. by the effective intensity of the immune response. 

The spatial pattern emerges when the environmental noise intensity a is 
properly tuned. Combined with the multiplicative noise, it should allow the 
system to perform transitions between two "phases". Too strong additive 
noise would however dominate the picture, preventing the formation of a 
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Fig. 10. Time-evolution of x, y, z and p in the system with additive Gaussian noise 
and diffusion, recorded at point [20, 20] on the spatial lattice. K\ = k\ — A = 
1.25, A = 0.5, D = 0.05, a = 0.01, A = 0, fc 2 = 0.1, fc 3 = 0.1. Initial conditions: 
x = 0,y = 0.4,z = 0,p = 0. 

pattern. Additionally, the diffusion parameter D should be carefully cho- 
sen, so that the pattern could have sufficiently distinct boundaries and be 
relatively stable. If the diffusion rate is high, the pattern dissolves quickly, 
whereas at small diffusion rate it only forms small "grains". The quantitative 
analysis of the interplay between the above-mentioned factors in the process 
of pattern formation merits a further study. 

After a sufficiently long time, the immune cells prevail globally. This 
turns out to be the effect of reflecting boundary conditions, which pre- 
vent the population densities from exceeding or 1. The existence of such 
boundaries causes that the trajectories of the system prefer to move in the 
direction of greater population of immune cells. A replacement of the ad- 
ditive Gaussian noise £(7*,^ with a multiplicative Gaussian noise, and a 
comparison with the model described here, would be another interesting 
issue for future research. 

The author would like to thank Dr. Ewa Gudowska-Nowak for helpful 
discussions and comments. 
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Fig. 11. An example trajectory of the system with additive Gaussian noise and 
diffusion, recorded at point [20,20] on the spatial lattice. The immune response 
parameter K\ = k\ + A = 2.25. The trajectory remains in the neighbourhood of 
the "y-phase" attractor and climbs towards maximal values of y. X = 0.5, 1? = 
0.05,cr = 0.01, A = 0, k 2 = 0.1, k 3 = 0.1. Initial conditions: x = 0,y = 0.4, z = 
0,p = everywhere. 
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Fig. 12. The time-evolution of x,y, z and p. in the system with additive Gaussian 
noise and diffusion, recorded at point [20, 20] on the spatial lattice. K\ = k\ + A = 
2.25, A = 0.5, D = 0.05, a = 0.01, A = 0,fc 2 = 0.1, k 3 = 0.1. Initial conditions as 
above. 
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Fig. 13. Trajectory of the system (@J recorded at point [20, 20] on the spatial 
lattice. Simulation time T = 4000. The trajectory jumps between two possible 
phases. Initial conditions: x — 0, y — 0.4,2 = 0,p = everywhere. Parameters: 
A = 0.5, D = 0.05, a = 0.01, A = 0.5, fa = 1.75, fa = 0.1, fa = 0.1, 7 = 0.01 




Fig. 14. Trajectory of the system |@J| recorded at point [20, 20] on the spatial 
lattice, y-z-projection, simulation time T — 10000. The trajectory jumps between 
two stable branches of stationary points, but, finally, it climbs up the "y-phase" 
branch. Initial conditions: x = 0, y = 0.4, z = 0,p = everywhere. Parameters: 
A = 0.5, D = 0.05, a = 0.01, A = 0.5, fa = 1.75, fa = 0.1, fa = 0.1,7 = 0.01. 
Curves show how the shape of branches of stationary points changes with switching 
value of fa ± A. Sloped line: the initial plane E — 0.4. 
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Fig. 15. Time evolution of x, y, z and p in system ijljl. Simulation time T = 10000. 
Initial conditions: x = 0,y = 0.4, z = 0,p = everywhere. Parameters: A = 
0.5, D = 0.05, a = 0.01, A = 0.5, fa = 1.75, fa = 0.1, fa = 0.1, 7 = 0.01 
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